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Abstract 

We provide explicit expressions for quadrature rules on the space of C 1 quintic 
splines with uniform knot sequences over finite domains. The quadrature nodes 
and weights are derived via an explicit recursion that avoids an intervention of 
any numerical solver and the rule is optimal, that is, it requires the minimal 
number of nodes. For each of n subintervals, generically, only two nodes are 
required which reduces the evaluation cost by 2/3 when compared to the clas¬ 
sical Gaussian quadrature for polynomials. Numerical experiments show fast 
convergence, as n grows, to the “two-third” quadrature rule of Hughes et al. 
m for infinite domains. 
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1. Introduction 

Numerical quadrature has been of interest for decades due to its wide ap¬ 
plicability in many fields spanning collocation methods [26], integral equations 
[3], finite elements methods [27] and most recently, isogeometric analysis [5]. 
It is also a preferable tool for high-speed solution frameworks din] as it is 
computationally very cheap and robust when compared to classical integration 
methods 1151 . 

A quadrature rule, or shortly a quadrature, is said to be an m-point rule , if 
m evaluations of a function / are needed to approximate its weighted integral 
over an interval [a, b\ 


oj(x)f(x) dx 




(1) 
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where ui is a fixed non-negative weight function defined over [a, h\. Typically, the 
rule is required to be exact , that is, R m .(f) = 0 for each element of a predefined 
linear function space C. Moreover, the rule is said to be optimal if to is the 
minimal number of nodes r, at which / has to be evaluated. 

In the literature, the term optimality may also refer to the approximation 
error that the quadrature rule produces. That is, the number of nodes is given 
and their layout is sought such that it minimizes the error for a given class of 
functions. Kohler and Nikolov H2GS] showed that the Gauss-type quadrature 
formulae associated with spaces of spline functions with equidistant knots are 
asymptotically optimal in non-periodic Sobolev classes. This is a motivation 
for studying Gauss-type quadrature formulae for spaces of spline functions, in 
particular, with equidistant knots. In this paper, by optimal we exclusively 
mean quadrature rules with the minimal number of nodes. 

In the case when C is the linear space of polynomials of degree at most 
2 to. — 1, then the TO-point Gaussian quadrature rule [T51 is both exact and opti¬ 
mal. The Gaussian nodes are the roots of the orthogonal polynomial n m where 
(•7To, 7Ti,..., 7r m ,...) is the sequence of orthogonal polynomials with respect to 
the measure n(x) = u>{x)dx. Typically, the nodes of the Gaussian quadrature 
rule have to be computed numerically, which becomes expensive, especially for 
high-degree polynomials. 

Naturally, getting more degrees of freedom is not achieved by using a higher 
polynomial degree, but rather by using polynomial pieces, i.e., splines. Addi¬ 
tionally to the polynomial case, the interval of interest is provided by a non¬ 
decreasing sequence of points known as a knot sequence (or vector), points where 
the resulting spline is considered to have a lower continuity. The knot sequence 
defines the local support of each basis function, that is, each one acts only lo¬ 
cally on a particular subinterval of the domain, and-depending on the knots 
multiplicities- spans a particular number of knots. We refer to [HI H2J [H] for a 
detailed introduction to splines. 

Regarding the quadrature rules for splines, Micchelli and Pinkus |19j derived 
the optimal number of quadrature nodes and specified the range of intervals, the 
knot sequence subintervals, that contain at least one node. There are two main 
difficulties compared to the polynomial case: firstly, the optimal quadrature rule 
is not in general guaranteed to be unique, e.g., when the boundary constraints 
are involved, and, secondly, jT9l determines only a range of intervals, i.e., each 
node has several potential subintervals to lie in. The latter issue is crucial 
as one cannot apply even expensive numerical solvers, because the algebraic 
system to solve is not known. For each assumed layout of nodes, one would 
have to solve a particular algebraic system using e.g., [1113 [55]. However, the 
number of eventual systems grows exponentially in the number of subintervals 
and therefore such an approach is not feasible. Instead, theorems that derive 
exact layouts of the nodes are essential. Our work in this paper contributes such 
a theoretical result for a particular space of splines. 

The quadrature rules for splines differ depending on the particular space of 
interest Sd, c , where d is the degree and c refers to continuity. For cases with 
lower continuity, the “interaction” between polynomial pieces is lower and hence, 
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naturally, a higher number of nodes is required for the optimal quadrature rule. 
The choice of the domain can bring a significant simplification. Whilst there are 
few rules that are exact and optimal over a finite domain, their counterparts are 
known when the integration domain is the whole real line [16] . The half-point 
rules of Hughes et al. are independent of the polynomial degree and the “half” 
indicates that the number of quadrature points is roughly half the number of 
basis functions. 

The half-point rules can be altered even for spaces with lower continuity, 
e.g., for £ 4,1 an optimal rule was also derived in [16]. The rule is called a “two- 
third rule” as it requires only two evaluations per subinterval whilst the classical 
Gaussian rule for polynomials needs three nodes. However, these rules are exact 
only over the real line. Because in most applications a finite domain is needed, 
additional nodes have to be added to satisfy the boundary constraints and a 
numerical solver has to be employed [4]. We emphasize that we focus here only 
on optimal rules as there are many schemes that introduce redundant nodes in 
order to overcome the problem with a finite interval, see [4] and the references 
cited therein. 

Regarding optimal rules over finite domains, Nikolov }22| proved the unique 
layout of nodes of the quadrature rule for 63,1 with uniform knot sequences and 
derived a recursive algorithm that computes the nodes and weights in a closed 
form. Recently pQ, the result was generalized for 5 3: 1 over a special class of 
non-uniform knot sequences, called symmetrically-stretched. The rules possess 
the three desired properties, i.e., they are exact, optimal and act in a closed 
form fashion, without intervention of any numerical solver. 

We emphasize that the computation of the nodes and weights of the optimal 
spline quadrature, also called Gaussian , is rather a challenging problem as one 
has to first derive the correct layout of the nodes and then, typically, to solve 
non-linear systems of algebraic equations. For higher degrees, the use of a 
numerical solver seems to be unavoidable. 

In this work, we derive a quadrature rule for spaces of C 1 quintic splines, S'sp, 
with uniform knot sequences over finite domains. The rule is exact, optimal, 
and-even though the degree is five-explicit. We also show numerically, when the 
number of subintervals grows, that the rule rapidly converges to the “two-third” 
rule of Hughes m- 

The rest of the paper is organized as follows. In Section [2] we recall some 
basic properties of and derive their Gaussian quadrature rules. In Section[3j 
the error estimates are given and Section [4] shows some numerical experiments 
that validate the theory proposed in this work. Finally, possible extensions of 
our method are discussed in Section [5] 

2. Gaussian quadrature formulae for C 1 quintic splines 

In this section we state a few basic properties of £ 5,1 splines and derive 
explicit formulae for computing quadrature nodes and weights for spline spaces 
with uniform knot sequences over a finite domain. Throughout the paper, 7 
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Figure 1: Four consecutive knots #/ c _ 2 ,..., of a uniform knot sequence, each of multiplic¬ 
ity four. Six non-normalized spline basis functions -D 4 fc— 3 ,..., -D 4&+2 with non-zero support 
on [xk— i,ajfc] are displayed. 


denotes the linear space of polynomials of degree at most d and [a, b] is a non¬ 
trivial real compact interval. 

2.1. C 1 quintic splines with unifoi'm knot sequences 

We detail several properties of spline basis functions. We consider a uni¬ 
form partition X n = (a = Xq, x±, ..., x n -i, x n = b) of the interval [a, b] with n 
subintervals and define h := A = Xk — Xk-\ for all k = 1 ,n. We denote 
by Slf , the linear space of C 1 quintic splines over a uniform knot sequence 
X n = (a = x 0 ,X!,..., x n = b) 

SI r = {/ € C^b] : f\ {xk _ 1>Xk) G 7r 5 ,fc = 1,..., n}. (2) 

The dimension of the space S'”is 4n + 2. 

Remark 1. In the B-spline literature mmM’ the knot sequence is usually 
written with knots’ multiplicities. However, in the isogeometric analysis litera¬ 
ture, see e.g., JlHf, the knot vector is usually split into a vector carrying the 
partition of the interval and a vector containing continuity information (knot 
multiplicity). As in this paper the multiplicity is always four at every knot, we 
follow the latter notation and, throughout the paper, write X n without multiplic¬ 
ity, i.e., x k < x k +i, k = 0,...,n-l. 

Similarly to CQ E2 > we find it convenient to work with the non-normalized 
B-spline basis. To define the basis, we extend our knot sequence X n with two 
extra knots outside the interval [a, b] in a uniform fashion 

X-i = xq — h and x„+i = x n + h. (3) 

The choice of x_i and x n+ i allows us to simplify expressions in Section [2T2j but 
this setting does not affect the quadrature rule derived later in Theorem |2.1[ 
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We follow [TU] and denote by D = {A}i=i~ 2 the basis of where 

P^Aksif) [*£fc — 25 X k — 25 X k — 1; X k — 11 Xk— 1 7 1 , Xk\ (■ 

P 4 k — 2 (t) — [^fc — 2, X k — 1, Xk — It Xk — 1, Xk— It X k , Xk\{- t)^_ , 

-^4fc —1 (t) — [^ft —1 t X k — 17 *Efc — 17 *£fc — 1 7 *^fc 7 7 ^fc] (‘ t) , 

P 4 k (t) — —1 7 ^fc — 15 *£fc — 15 X k , Xk 7 *£fc 7 (• t)-j- • 


where [.]/ stands for the divided difference and it + = max(it, 0) is the truncated 
power function. The direct computation of the divided differences gives the 
following explicit expressions for t £ [xk- 2 , x k -\\ 


-D4/C-3W 

D 4 k-2{t) 


and for t £ [xk-i,Xk] 


D^ksit) 
P 4 k—2 if) 
D 4k -i(t) 
D 4 k{t) 


(t-xfc_ 2 ) 4 (ccfc+8a:fc_i-9t) 

4ft 6 

(t-Xk-j) 5 

4ft 6 > 


(x k -tf 

4ft 6 ’ 

(z fc -Q 4 (:r fc _2+8:r t ,_i-9t) 
4ft 6 ’ 

10(t-x fc _i) 2 (a:fc-i) 3 
ft 6 > 

10(t-Xfc_i) 3 (ccfc-i) 2 

ft 6 


(5) 


( 6 ) 


The functions have the following pattern: six basis functions D 4 k- 3 , ..., £* 4^+2 
have non-zero support on [£*,_!, a:*,], moreover, two of them, D 4k -\ and D 4k , 
act only [ Xk-i,Xk\ and are scaled Bernstein basis functions, see 0 and Fig. [I] 
Among the basic properties of the basis D, we need to recall the fact that 
D 4 k+ 2 {t) < D 4 k+i{t) for t € [ Xk-i,Xk] for k = 1 and are equal at the 

knot x k , that is, D 4 k+i(xk) = D 4k + 2 {xk) = jr- Moreover, we have that 
D 4k -i ( Xk ~ 1 2 +Xk ) = D 4k ( Xk ~ 1 2 +Xk ) = From 0 and 0, the integrals of 

the basis functions are computed 


I[Dk\ = g for k = 3,4, ..., An, 


(7) 


where /[/] stands for the integral of / over the interval [a, b\. The first and the 
last two integrals are 

I[D 1 ]=I[D An+2 ] = ± and I[D 2 ] = I[D 4n+1 ] = i. (8) 

A less obvious property that binds together four consecutive basis functions, 
which is used later for our quadrature rule in Section [272] is formalized as follows 

Lemma 2.1. Let X n = (a = Xq,x 4 , ...,x n = b), be a uniform knot sequence and 
for any k = 1,..., n define 

P k [t) = 2L> 4 ft+i(t) — 2D 4k +2(t) + -D 4 k-\{t) — D 4 k(t). (9) 


Then Pk(t) > 0 for any t, S (xk-i,Xk) and P k {t) = 0 if and only ift= Xk ]^ Xk 
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2£>4fc+l ~ 2I?4fc + 2 



Figure 2: The linear blend of basis functions & Pfc, is expressed as a Bezier curve on 

& _ 1 , Xk ~ 1 + Xk ] w ith a control point sequence (red dots) with non-negative ^/-coordinates 
. Consequently, P ^ is non-negative on [xk_i,Xk\ and has a single root (of multiplicity 
two) at 


Proof. Over an interval ( Xk-i,Xk ), the function Pk is a single quintic polyno¬ 
mial. Therefore, it can be expressed in terms of Bernstein basis and can be 
viewed as a Bezier curve on a particular domain, see Fig. [2] Looking at its 
shape, one cannot conclude non-negativity from the control polygon, when con¬ 
sidered on the whole interval ( Xk-i,Xk )• Hence we define 


Pk(t) = Pk(t ) on [sfc-i, Xk ^ Xk 
p k(t) =Pk(t) on [ Xk ~ 1 2 +Xk ,x k ], 


( 10 ) 


and using h = Xk — %k-i we further write 


P^(t) = ZtfPm, where B?(t) = (?) (^^^ ( 2 *»-*+*- 2 « ) 


5 -* 


and analogously for P%. The conversion from monomial to Bernstein basis gives 
the control points (pj,... ,p\) of P£ over the interval [xk-i,Xk~i + h/2 ] as 


(P0.P1.P2.P3.P4.Ps) = ^ 0 . 0 , ^, 0 , 0 ^ 

and similarly for P| we obtain 

(p^p?,pI,p^p!,^) = (o,o, i ^,^,^,o) 


( 11 ) 


( 12 ) 


Therefore, P^ and P^ are non-negative on open intervals (xk-i,Xk-i+h/2) and 
(. Xk-i + h/2,Xk), respectively. Due to the fact that (p\,p\) = (Po.Pi) = (0.0). 
the only root (with multiplicity two) of Pfc on ( Xk-i,Xk ) is Xk-i + h/2. □ 
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ZD^k+i — 2D 4fe+2 



Figure 3: The assumption of ex iste nce of a single node r % inside [xk,%k+l] implies r, = 
^fc+i+^fc Consequently, the rule |l3| must return zero for Pj. on [x/._ ], X/,.] . i.e. Q%^_ 1 {Pk) = 

0. As P/ :: is non-negative on (xfe_ \ ,Xk) with one double root Xk ~ 1 + Xk this fact violates the 
assumption of a single node in (xk,Xk+\). 


Remark 2. D^k+i —D^+ 2 , D^k-i, D^k o.re all positive polynomials on (xk-i, Xk) 
and therefore there exists infinitely many non-negative blends. However, the ex¬ 
istence of a non-negative blend when the coefficients have to satisfy a certain 
constraint is not obvious and the full impact of this particular blend with coeffi¬ 


cients 2 


|,-i 


will be seen later in Lemma 


2 . 2 . 


2.2. Gaussian quadrature formulae 

In this section, we derive a quadrature rule for the family Sf j. see (|2j). 
Similarly to |Tj, we derive a quadrature rule that is optimal, exact and explicit. 

With respect to exactness and optimality, according to mm there exists 
a quadrature rule 

i>b 2n+l 

Qa(f) = / (13) 

Ja i=1 

that is exact for every function / from the space Sf t . The explicitness follows 
from the construction. 

Lemma 2.2. Let X n = (a = Xq, Xi, ..., x n = b) be a uniform knot sequence. 
Each of the intervals Jk = (xk-i,Xk) (k = 1,... ,n) contains at least two nodes 
of the Gaussian quadrature rule \13 1 ). 

Proof. We proceed by induction on the index of the segment Jk- There must 
be at least two nodes of the Gaussian quadrature rule inside the interval J±. If 
there were no node inside Ji, the exactness of the rule would be violated for Di. 
If there was only one node, using the exactness of the quadrature rule for 
and -D 4 , it must have been the midpoint t\ = Go+xi) w j^j 1 the weight oq = j^h. 
However, this contradicts exactness of D\ and D 2 as D\ < D 2 on {xq,Xi). 

Now, let us assume that every segment Jk, k < n, contains-two or more- 
Gaussian nodes and prove that Jk+i contains at least two nodes too. By con¬ 
tradiction, if there is no node inside (xk, Xk+ 1 ), the exactness of the quadrature 
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oik 


Xk— 1 


-o — 

T2k-1 


T~2k 


'X k 


Figure 4: Notation on [xj._i,a;*.]. 


rule (13) for D A k +3 is violated. If there is a single node in ( Xk,Xk+i ), due to 


the exactness of the quadrature rule (13) for £> 4^+3 and D±k+i, it must be the 
midpoint = P k+Xk +J as it j s their only intersection point, see Fig. [ 3 J and 
their integrals are equal, see (|7|, I[D Ak + 3 \ = I[Dik +4] = =. 


ing weight must be Wj = ^ h. 
combining with (Vi, we have 


The correspond- 


Moreover Z> 4 fc+i(^) - D ik + 2 {n) = - 


el h and 


2wi(D 4 fe+i(n) - D Ak +2{xi)) = - 


12 ' 


(14) 


Consider a blend Pk{t) = 2D A k+i(t) — 2,D A k+2{t) + A D A k-i{t) — D A k(t), see 


Fig. i As Pk is a blend of basis functions, the rule (jT3j) must integrate it 
exactly on [xk-i, Xk+i], that is Qxtt\(Pk) = I(Pk) = — pj. However, combining 


this fact with (14), the rule must return zero when applied to Pk on [xk-i, Xk\ : 

Pk is non-negative on ( Xk-i,Xk ) 


i.e. 


Q% k k _ 1 (,Pk) = 0. But due to Lemma 


with the only root at Xk ]^ Xk , which contradicts the assumption of a single 
quadrature node in (xk,Xk+ 1 ) and completes the proof. □ 


2.1 


Corollary 1. If n is an even integer , then each of the intervals Jk = ( Xk~i,Xk ) 
(fc = 1,2,..., n) contains exactly two Gaussian nodes and the middle x n / 2 = 
(a + b)/2 of the interval [a,b\ is also a Gaussian node. If n is odd then each 
of the intervals Jk = (xk-i,Xk) (k = 1,2,... ,n;k ^ (n + l)/ 2 ) contains exactly 
two Gaussian nodes, while the interval J( ra +i )/2 contains three Gaussian nodes: 
the middle (a + b)/2 and the other two positioned symmetrically with respect to 
(a + h)/2. 


Proof. From [15] , 
Gaussian nodes. 


From Lemma 2.2 


the optimal quadrature rule (131 is known to require 2n + 1 


we know the location of 2 n of them as 
each of the intervals Jk contains at least two nodes. The last node must be 
the midpoint ( a + b)/2. We prove it by contradiction, distinguishing two cases 
depending on the parity of n. For n even, if one of the intervals Jk has more 
than two nodes then, by symmetry, J n -k has to contain the same number of 
nodes and we exceed 2n + 1, contradicting our quadrature rule (13). For n 


odd, let us assume that the middle interval J(„+i )/2 contains exactly two nodes. 
Then, by symmetry, at least two of the remaining intervals contain three nodes, 


contradicting our quadrature rule (13). Therefore, the middle interval J( n + 1)/2 
contains exactly three nodes, where the middle one is, again by symmetry, forced 
to be the midpoint (a + b)/2. □ 


With the knowledge of the exact layout of the optimal quadrature nodes, we 
now construct a scheme that starts at the boundary of the interval and parses to 


















its middle, recursively computing the nodes and weights. This process requires 
to solve only for the roots of a quadratic polynomial. Let us denote 

Oik — X 2k — 1 X k — 1) 0k — Xk X 2 fc; (lb) 

where r 2k -i and r 2k , T 2k ~1 < T 2 fc, are the two quadrature nodes on ( Xk-i,Xk ), 
k = 1,..., [n/2] + 1, see Fig. |4j Keeping in mind h = Xk — x k ~ i, we have 

Xk T 2 k—1 h Ofc, X 2 k Xk— 1 h 0k * (16) 


Let w 2 fc-i and W 2 fc be the corresponding weights of the Gaussian quadrature rule 
over the interval (x k -i,Xk)- The exactness requirement of the rule when applied 
to D 4 k-i and D Ak , see § and ([7]), gives the following algebraic constraints 


uj 2k -ia 2 k (h - a k ) 3 + u 2k (h - 0 k ) 2 01 = ^, 

h 6 

u] 2k -ia 3 k (h - a k ) 2 + u} 2k {h - 0 k ) 0l = xx- 


(17) 


The exactness of the rule when applied on and D Ak - 2 , respectively, gives 


W2fc—1 


/ 5 (h - a fc ) 4 

V 2 h 5 


u 2 k-i(h 


9{h - a fc ) 5 \ 
4 h 6 ) 


+ X> 2 k 


o: k ) 5 + u 2k 0k 

, 2 h 5 4 h 6 ) 


^h?r 4fc-3) 


^4fc-2, 


(18) 


where r 4*,_3 and ?' 4 fc _2 are the residues between the exact integrals, see 0 and 
Q, and the result of the rule when applied to D Ak -3 and D^k -2 on the previous 
interval [ 2 ^- 2 , x k -i\, respectively. That is 


XAk -3 = 7 [D 4fe _ 3 ] - QZ k k _\{Dik-3), 
XAk-2 = I[DAk-2\ ~ Q-xtll i^Ak-2)- 


(19) 


Due to the fact that both ( |17| and (18) are linear in u) 2k -i and w 2 fc> their 
elimination from gives 


X>2k -1 — 


W2 k — 


h 5 (h - 20 k ) 


60 a 2 k (h - a k ) 2 (h - a k - 0k) 1 
h 5 (h — 2a k ) 

60/3 l{h - 0 k ) 2 {h -a k - 0k) ’ 


and from (18) we obtain 

-2h 5 (90 k r4k-3 - 10/ir 4 fe- 3 + 0kX4k- 2 ) 


W2fc-1 — 


W2fc = 


5 (h - a k ) 4 (h - a k - 0 k ) 

—2h 5 (hr4 k -3 + a k r4k-2 + 9r 4 fc-3«fc ~ hr Ak _ 2 ) 
50%(h-a k -0k) 


( 20 ) 


( 21 ) 
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Figure 5: The algebraic system ( |22| > over the domain (0, h) x (0, h) (grey) for the first (k = 1) 
subinterval [a;o,®i] is shown. The two intersection points correspond to the two Gaussian 
nodes on [xo,xi\ and are computed by projection onto ai-axis using the resultant. The 
coordinates of the intersection points with respect to a ±-axis are the roots of the quadratic 
polynomial \2A\ . 


Equating 0 J 2 k-i 


and uj 2 k from (20) and we obtain 


®k(ak,/3k) = 0 , 

^k{otk,Pk) = 0 ; 


( 22 ) 


an algebraic system of degree three with the unknowns a k and (3 k ■ Solving this 
non-linear system of two equations with two unknowns can be interpreted as the 
intersection problem of two algebraic curves, see Fig. [5} The domain of interest 
is (0 ,h) x (0 ,h) as both quadrature points lie inside (xk-i,Xk)- 

Using the resultant, see e.g. [9), of these two algebraic curves in the direction 
of /3k, one obtains a univariate polynomial, in general, of degree nine. Interest¬ 
ingly, our system ( |22| produces-for all admissible residues and r 4 fc_ 2 -only 

a quintic polynomial Ek{a k ) that gets factorized over R as 

^CS , f3k), ^ki^k, Pk ) ) -Ufc ) Qk (ckfc ) , (23) 


where Qk is a quadratic factor and the vector of its coefficients with respect to 
the monomial basis, cft no = (<?o, , q 2 )> is-in terms of the residues-expressed as 

92 = 1 - 480r 4 fc_ 3 + 576rl k _ 3 + 7>7Qr 2 ik _ 2 - 1152r 4 fc- 2 r 4fc _3, 

q* = 2h{l2r Ak _ 2 + 108r 4fc _ 3 - 1), (24) 

9o = ft2 (! ~ 24r 4 fe-2 + 24r 4fc _ 3 ), 

and for the vector of monomial coefficients = (cq, c\, c 2 , c 3 ) of the cubic 
factor C k we obtain 


c 3 = —216r4fc_3 — 24r4fc_2 + 2, 

c 2 = /i(24r 4 fe_2 - 24r 4 fc_ 3 - 5), 
c\ = Ah 2 , 
c o = ~h 3 - 


(25) 
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%m —2 

-ho- 


%m— 1 

4 - 


x n 

-th 


^m+1 


■o+ 


7n+l ^”n+2 


Figure 6: The situation for odd n on the middle interval J m = [xr, 
is the middle of the interval, r n and r n _|_2 are computed from (|32|). 


-l ,Xm\- The node r n +i 


A Maple worksheet with this algebraic factorization is attached to this submis¬ 
sion. 


We recall that two roots of E k (23) determine the two quadrature nodes that 


lie inside [x k -i,Xk\. Interestingly, the cubic factor does not contribute to the 
computation of the nodes which is formalized as follows. 

Lemma 2.3. The cubic polynomial Ck p?5| ) has no roots inside [0, h]. 

A proof can be found in Appendix. 

We now proceed to the main contribution of the paper, a recursive algorithm 
that computes the nodes and the we ight s of the Gaussian quadrature for uniform 
C 1 quintic splines. Due to Lemma 2.3 the recursion operates in a closed form 


fashion by solving for the roots of quadratic polynomial (24). Before we state 
the theorem, we need to establish some notation. 

Let us denote by A k and B k . k = 2,..., [ra/2] + 1, the actual values of 
residues (191 when being evaluated at the nodes T 2 k -3 and r 2 fc _2 on the interval 
[xk- 2 , Xk-i], he., 


Ak = I[D^k-3] — W2fc-3-D4A:—3(T2/c—3) _ a/2fc-2-D4/c-3( 72 / 0 - 2 ), 
Bk = h[Z?4 fe _ 2 ] — W 2 fc-3^4fc-2 (t2/c— 3 ) — W 2 fc_2D4 fc _ 2 (r 2 fe_ 2 ), 


(26) 


and the coefficients of the quadratic polynomial (241 become 


(27) 


a k = 1 — 480A fc + 576A^ + 576 Bl - 1152 B k A k , 
b k = 2h(YlB k + 108A fc - 1), 
c/o = h~ (1 24Bfc + 24 A k )- 

In the case when n is odd, see Fig. [6j the middle subinterval contains three 
nodes and the algebraic system that needs to be solved results in a quadratic 
polynomial with the following coefficients 

a m = —2(108A m + 12B m + 1), 

b m = 2/i(108A m + 12 B m - 1), (28) 

Cm = h 2 (12A m — 12 B m + 1). 

We are now ready to formalize the main theorem. 

Theorem 2.1. The sequences of nodes and weights of the Gaussian quadrature 


rule (13) are given explicitly by the initial values A k = A an d B k = | and the 


recurrence relations (k = 1,..., [n/2]) for the nodes 


72/0-1 — %k-l 


-bk-y/bl~4:a k c k 

2a k 


and T 2k = x k — 


-b k + ^/bl-4a k c k 


2 a/e 


(29) 
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and for the weights 


_ _ h 5 (2r 2k - 2x k -i - h) _ 

L ° 2k ~ 1 60(r 2fc _i - x k + h) 2 (x k - T 2 fc-i) 2 (r 2 fc - T 2 k -i ) ’ 

_ _ h 5 (2x k - 2r 2fc _i - h) _ 

60(a: fc _i - r 2k + h) 2 {T 2k - a; fe _i) 2 (r 2fc - r 2fe _i)' 

If n is even (n = 2m) then r n +i = x m = (a + 6)/2 and 


(30) 


w n +i — 4 h( - — A m+ i). 
6 

If n is odd, (n = 2m — 1) then r„+i = (a + b)/2, 

— b m — 'fbi^— a ta rn 'c 


T~n — %m—l T 


and r n+2 = x m - 


(31) 


—bm+y/b^—Aa^Cm (32) 


and the corresponding weights are 

_ h(108A m + 12 B m - l) 2 
x n — u] n+2 — 


XI n +1 — 


30(156A m - 36B m + 1) ’ 

Ah{1152A m B m + 264 A m - 576 A 2 m - 57 QB 2 m - 24 B m + 1 ) 
15(156A m - 36+ 1) 


(33) 


Proof. We proceed by induction. Assume the quadrature nodes [t 2 i_i,t 2 i) and 
weights (oj 2 i-i,uj 2 i) are known for l = 1,..., k — 1 (k < [n/2]) and compute the 
new ones on (x k -i,x k ). For k = 1, as there are no nodes on [x-\,xq), (26) 
gives A\ = I[Di] = and B\ = I[D 2 \ = |. By Corollary [I] there are exactly 
two nodes inside (x k -i,x k ). Due to Lemma 2.3 only the roots of the quadratic 
factor in (23) contribute to the computation of the nodes and hence solving 


Qk{oi k ) = 0 with coefficients from (27) gives a k and f3 k . Combining these with 


(15) results in ([29]). The weights are computed from (20) using the identities 
(15) and ( |l6| ). By Corollary [lj we have r n +i = (a+ b)/2. If n is even, using the 
exactness of the quadrature for D 2n+ i, the associated weight is computed from 


g — I[D 2n + i] — A m+ 1 + ui n +iD 2n+ i(T n+ i). 


(34) 


Evaluating D 2n+ i((a + &)/2) = ^ gives (31). If n is odd, due Corollary [lj there 
are three nodes inside (a; m _i,a; m ); one is the middle point (a + b)/2 and the 
other two are symmetric with respect to it, see Fig. [6] Using the notation of 
(15) for the middle interval, i.e., a m = r 2m _i — x m _i, the rule must integrate 


exactly D± m _ 2 and D^ m _i which gives the following 3x3 algebraic 

system 

(h— a m ) 5 +a^, , . , 1 _ ^ 

Ih B UJn ' 128 h U>n + 1 


( h—ctm ) (9a m +/i)+a(^ (lOh—9a m ) 


i 11 . , 

"t" 128/i Xn+l 


4 


= B n 


10 a((/t-a n 

W 


4- _JL, , — I 

' ^e>h u 'n+l 


(35) 
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Algorithm 1 


GaussianQuadrature([a, 6], n ) 


1: INPUT: compact interval [a, b] and number of uniform segments n 


2 

3 

4 

5 

6 

7 

8 
9 

10 

11 

12 

13 

14 

15 

16 


Ai 24 ; I?i — g ; 
for k = 1 to [n/2] do 

compute r 2k _ i, r 2k from ([29]), and u 2k -i, u 2k from ( |30| ); 
end for 

T n+ i = (a + b)/ 2; /* middle node */ 

if n is even then 

compute w n +i from (31); 
else 


compute T n and r n + 2 from (32) and oo n +i and oj n + 2 from (33); 

end if 

for k = 1 to [n/2] do 


T~2n-2k+3 — T 2k -1] T 2n -2k+ 2 — T 2k ] 


/^symmetry */ 


^>2n-2k+3 — W 2 fc-i; ^>2n-2k+2 — <^2k] 

end for 

OUTPUT: {r i5 set of nodes and weights of the Gaussian quadra¬ 

ture on interval [a, &]; 


with unknowns a m , uj n and w n +i. Eliminating from the first two and 

second two equations, respectively, and solving for ui n we obtain 

= 2 h 5 (llA m -B m ) = /i 5 (240f? m — 11) 

5(/i 2 - 2ha m + 2a‘? n )(h - 2a m ) 2 60(h 2 + 9ha m - 9a^)(/i — 2a m ) 2 


and the problem reduces to solving for the roots of a univariate (rational) func¬ 
tion in a m . The numerator is a quadratic polynomial with coefficients (28) 
which proves (32). Inserting (28) into (35) and solving for oj n and w n +i then 
gives (33) and completes the proof. □ 

For the convenience, we summarize the recursion in Algorithm [l] 


3. Error estimation for the C 1 quintic splines quadrature rule 

In the previous section, we have derived a quadrature rule that exactly inte¬ 
grates functions from 1 with uniform knot sequences. If / is not an element 
ofSj?!, the rule produces a certain error, also known as remainder. The analysis 
of this error is the objective of this section. 

Let W[ = {/ £ C r ~ 1 [a,b\- 1 /( r_1 )abs. cont., H/^Hjq < oo}. As the 
quadrature rule © is exact for polynomials of degree at most five, for any 
element / £ Wf. d > 6, we have 

R2n + l[f\ ~ I[f] - Q b alf ] = [ b K 6 (R 2n+1 ;t)f^\t)dt, 

J a 
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where the Peano kernel m is given by 


Ke{R 2 n+ljt) — i?2n+l 


5! 


An explicit representation for the Peano kernel over the interval [a, b] in terms 
of the weights and nodes of the quadrature rule (131 is given by 


2n+l 


K 6 (R 2 n+i-,t) = -^ u k {t-T k )\. 


720 


(36) 


k =1 


Moreover, according to a general result for monosplines and quadrature rules 
m , the only zeros of the Peano kernel over (a, b) are the knots of multiplicity 
four of the quintic spline. Therefore, for any t G ( a,b ), I\Q(R, 2 n +i', f) > 0 and, 
by the mean value theorem for integration, there exists a real number £ G [a, b] 
such that for / € C 6 [a, b} 

[■b 

i?2n+l(/)=C 2 ra+1,6 / (6) (0 with Q2n+1,6 — / -^6(-^2n+lj tydt. (37) 

J a 

Hence, the constant C 2 n +i ,6 of the remainder i? 2 n+i is always positive and our 
quadrature rule belongs to the family of positive definite quadratures of order 


6, e.g., see mmm- Integration of ([36]) gives 


Theorem 3.1. The error constant C 2 n +i ,6 *n (5?) of the quadrature rule (13) 
is given by 


2n+l 


_ (6 — a ) 7 1 ^ . , 6 

C2n+1,6 — cr ,^ n 79n 2^ a ) ■ 


5040 


(38) 


fc=i 


4. Numerical Experiments 

In this section, we show some examples of quadrature nodes and weights for 
particular numbers of subintervals and discuss the asymptotic behaviour of the 
rule for n —> oo. 

In the case when the domain is the whole real line, the exact and optimal 
rule is easy to compute. Similarly to [16] . Eq.(29), where the rule was derived 
for SAi, for 5 5]1 case one obtains 


/(f) dt = 






(39) 


ie z 


that is, the nodes are the knots and the middles of the subintervals. Similarly 
to d, only two evaluations per subinterval are needed which gives 2/3 cost 
reduction ratio when compared to the classical Gaussian quadrature for poly¬ 


nomials. Observe the convergence of our general uniform rule to its limit, (39), 
when n —> oo. The weights and nodes are shown in Table [T] Only few initial 
nodes and weights differ from the limit values as ^ = 0.46 and ^ = 0.53. 
From Table |T] we conclude that for large values of n, one needs to compute only 
the first nine nodes and weights that differ from the limit values by more than 
e = KT 16 . 
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odd 


even 






Figure 7: The quadrature rules eu for various n are shown. The interval is set [a, b] = 
[0,n], i.e., the distance between the neighboring knots is normalized to h = 1. The green 
dots visualize the quadrature; the ^-coordinates are the nodes and the y-coordina 1 .es the 
corresponding weights. As n —> oo, the nodes converge to the knots and the midp oint s of the 
subintervals, the weights converge to 0.46 and 0.53 (black lines), cf. TablcJTjand |39|. 


5. Conclusion and future work 

We have presented a recursive algorithm that computes quadrature nodes 
and weights for spaces of quintic splines with uniform knot sequences over finite 
domains. The presented quadrature is explicit, that is, in every step of the 
recursion the new nodes and weights are computed in closed form, without using 
a numerical solver. The number of nodes per subinterval is two and hence the 
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Table 1: Nodes and weights for Gaussian quadrature ( |13|) with double-precision for 
various n are shown. To observe the convergence to ( |39| >, the interval was set as 
[a, b] = [0,n]. Due to the symmetry, only the first n + 1 nodes and weights are 
displayed. 



n = 

= 5 

n = 

= 6 

i 

Ti 

UJi 

n 

Ui 

i 

0.1225148226554413 

0.3020174288145723 

0.1225148226554413 

0.3020174288145723 

2 

0.5441518440112252 

0.4850196082224646 

0.5441518440112252 

0.4850196082224646 

3 

1.0064654716056596 

0.4467177201362911 

1.0064654716056596 

0.4467177201362911 

4 

1.5002730728687338 

0.3303872093804185 

1.5002730728687338 

0.5330387209380418 

5 

2.0000387957905171 

0.4665398664562177 

2.0000387972956304 

0.4665398713719121 

6 

2.5 

0.5333333108648244 

2.5000000105321137 

0.5333333220982075 

7 

- 

- 

3 

0.4666666568370204 


n = 

= 7 

n = 

= 8 

1 

0.1225148226554413 

0.3020174288145723 

0.1225148226554413 

0.3020174288145723 

2 

0.5441518440112252 

0.4850196082224646 

0.5441518440112252 

0.4850196082224646 

3 

1.0064654716056596 

0.4467177201362911 

1.0064654716056596 

0.4467177201362911 

4 

1.5002730728687338 

0.5330387209380418 

1.5002730728687338 

0.5330387209380418 

5 

2.0000387972956304 

0.4665398713719121 

2.0000387972956304 

0.4665398713719121 

6 

2.5000000105321137 

0.5333333220982075 

2.5000000105321137 

0.5333333220982075 

7 

3.3.00000000150452 

0.4666666617518435 

3.0000000015045293 

0.4666666617518435 

8 

3.5 

0.5333333333333333 

3.5000000000000000 

0.5333333333333333 

9 

- 

- 

4 

0.4666666666666665 


n = 

= 9 

n = 

10 

1 

0.1225148226554413 

0.3020174288145723 

0.1225148226554413 

0.3020174288145723 

2 

0.5441518440112252 

0.4850196082224646 

0.5441518440112252 

0.4850196082224646 

3 

1.0064654716056596 

0.4467177201362911 

1.0064654716056596 

0.4467177201362911 

4 

1.5002730728687338 

0.5330387209380418 

1.5002730728687338 

0.5330387209380418 

5 

2.0000387972956304 

0.4665398713719121 

2.0000387972956304 

0.4665398713719121 

6 

2.5000000105321137 

0.5333333220982075 

2.5000000105321137 

0.5333333220982075 

7 

3.0000000015045293 

0.4666666617518435 

3.0000000015045293 

0.4666666617518435 

8 

3.5000000000000000 

0.5333333333333333 

3.5000000000000000 

0.5333333333333333 

9 

4.0000000000000000 

0.4666666666666666 

4.0000000000000000 

0.4666666666666665 

10 

4.5 

0.5333333333333333 

4.5000000000000000 

0.5333333333333333 

11 

- 

- 

5.0000000000000000 

0.4666666666666666 


cost reduction compared to the classical Gaussian quadrature for polynomials is 
2/3. We have also shown numerically that in the limit, when the length of the 
interval goes to infinity, our rule converges to the “two-third rule” of Hughes et 
al. (16] that is known to be exact and optimal over the real line. 

For C 1 splines, our quadrature rule is optimal (Gaussian), that is, it re¬ 
quires minimal number of evaluations. However, it is still exact for any quintic 
splines with higher continuity and therefore we believe it can be used for various 
engineering applications. 


16 










































As a future work, optimal rules for spaces of higher polynomial degree and 
primarily of higher continuity are within the scope of our interest. Since the 
higher continuity reduces the number of optimal nodes, their layout in subin¬ 
tervals is difficult to assign and hence finding these is still an open problem. 
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Appendix 


Proof of Lemma \2.3[ Without loss of generality, we may assume h = 1 as the 
roots of Ck change with scaling factor h, cf. (25). The cubic Ck can be split 
into two summands fk and gt, the first independent and the latter dependent 
on r 4fc _ 3 and r 4fc _ 2 , he. 


f k = 2f 3 - 5f 2 + At - 1 = (t - 1 )\t - i), 

= (—216?' 4 fc_3 — 24r 4 fc_2)t 3 + (24r 4 fc_ 2 — 24r 4 fc_3)t 2 (40) 

= — 24f 2 (f - r 4fc-2 ~ ? 4k—3 j 
9r 4fc _ 3 + r 4 fc _ 2 


We show that Ck(t) < 0 on [0,1]. We denote by the non-zero root of gk, 


_ r 4 fc_ 2 — Tik -3 

9 9r 4 fc_ 3 + r 4 fe_ 2 ’ 

and consider the particular subdivision of [0,1] into 1) [0,£ g ], 2) [£ 9 , |], and 3) 
[i, 1]. We investigate Ck on each of these three subintervals separately: 

Case 3) We express Ck in Bernstein (B) basis and show all its coefficients 
are negative. Let Tji ^ be the transformation matrix m from monomial to 
Bernstein basis on [4,1] 


%!, - 


C 1 n\ 

i25i 
2 3 6 L 

1 _5_ 2 i 
4 12 3 

\I I I 1 / 

\ 8 4 2 / 


(41) 


and let c B = (cf , cf , cf , cf) be the vector of Bernstein coefficients. Then the 
conversion is given by c k B = c^ 0 Tji ^ and we obtain 


cf — — 33 r 4 fc _3 + 37 ' 4 fc_ 2 , 

c ? = Y 2 ~ 64r 4fc-3 + 4r 4fc _ 2 , 
cf = -124r 4fc _ 3 +4r 4fe _2, 
cf = —240r 4 fc_ 3 . 


(42) 


Looking at the cf coefficients, we start with the second one and prove that 
64r 4 fc_3 — 4r 4 fe _ 2 > fj. We know that r 4 fc _3 and r 4 fc _ 2 are by definition both 
positive and also r 4 fc _ 2 > r^ks, which is a direct consequence of D±k+i(t) > 
-D 4 fe+ 2 (f) on {xk~i,Xk), and also 7 ' 4 fc_ 3 ,r 4 fc _ 2 < |. Consider the blend Pk(t) 
from Lemma pLl) Pk(t) > 0 gives polynomial inequality 2D±k+i{t) — 2D 4 fc _|_ 2 (f)+ 
^D^k-iit) > Dik(t)- By evaluating both sides at the two nodes and by multi¬ 
plying with corresponding weights, i.e., by applying the quadrature rule Q on 
[ccfc_i,Xfc], we obtain 

2(4r 4 fc_3 — 7’ 4 fc- 2 ) + tw > v (43) 

12 0 
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because D^k-i and D 4 fc act only on this interval and hence the rule reproduces 
their integrals exactly. Combi ning (43) with r 4 fc _3 > 0 proves the desired in¬ 
equality. Moreover, combining (43) with | > r^k -2 gives 


16r 4fc _ 3 > 5r 4 fe- 2 , 


(44) 


and the other three inequalities follow directly from (44) and the fact r 4 ^_3 > 0. 

Case 1) Similarly to case 3), we compute the cf coefficients and show that 
all are negative. The conversion is given by c k B = cf io 7j 0i t ]. where 


T [o,fI 9 ] = 


(l l 1 1\ 

0 3 fg 5 £3 

0 i £ 2 £ 2 
0 0 


0 

Vo 


(45) 


is the corresponding transformation matrix. We obtain c k B = (cf, cf, cf, cf) 


= -l, 


C 3 = 


—31r 4 fe_ 3 + r 4 fc -2 
3(9r 4fe _ 3 +r 4fc _ 2 )’ 

4 80r 2 fe _ 3 - 5 r 4 fe _ 3 r 4fe _2 + 6 (r 4fe _ 3 - r 4fc _ 2 ) 3 
(9r 4fc _ 3 + r 4fe _ 2 ) 2 
— 100 r| fc _ 3 (llr 4 fc _ 3 - r 4 fe_ 2 ) 


(46) 


(9r 4fe _ 3 + r 4 fe _ 2 ) 3 


Negativity of cf and cf follows directly from (441. It remains to prove 
80r 4fe _ 3 - 5r 4 fe _ 3 r 4fc _ 2 + 6 (r - 4fc _ 3 - r 4fc _ 2 ) 3 > 0 , 


which using (441 and r 4 ^_2 > r 4 fc _3 simplifies to 


64r 


4fc-3 > ® r 4fe-2! 


which again follows from (44) using r 4 fc_ 2 < g. 

Case 2) fk has a double root at t = 1 and has a double root at t = 0. For 
the third root, it holds 


fik-2 — ^4fc-3 

9 r 4fe _3 + r 4fc _ 2 


< 


1 


(47) 


which follows from (44). From case 3), we know gk{\) = —33r 4 fc_3 + Sr^k -2 
and from case 1) fklfg) = cf that are both negative. Moreover, we know 
that fk is monotonically increasing while gk is monotonically decreasing on 
(£ g , ^). Therefore fk(t) < 0 on [£ g , |) and gk < 0 on (£ g , |] which completes 
the proof. □ 
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